Code
model {
# Specify priors
## Dynamics Parameters
for (t in 1:(n_years-1)){
phi[t] ~ dnorm(0, 0.1) #survival
gamma[t] ~ dnorm(0, 0.1) #colonization
}
## Occurrence-level Parameters
alpha01 ~ dnorm(0, 0.1) #Year 1 occurrence model intercept
for(n in 1:n_xcovs){
alphas[n] ~ dnorm(0, 0.1) #Occurrence covariate intercepts
}
## Detection-level Parameters
beta0 ~ dnorm(0,0.1) #beta0
beta1 ~ dnorm(0,0.1) #tmin covariate
beta2 ~ dnorm(0,0.1) #daylight covariate
beta3 ~ dnorm(0,0.1) #clut0 covariate
beta4 ~ dnorm(0,0.1) #Clut1 covariate
beta5 ~ dnorm(0,0.1) #Clut2 covariate
beta6 ~ dnorm(0,0.1) #Clut3 covariate (ref is 4)
beta7 ~ dnorm(0,0.1) #water indicator covariate
# Fit our Occurrence-level submodel
## First Year
logit_psi[1:n_sites, 1] <- alpha01 + xmat %*% alphas
## Following years
for (t in 1 : (n_years - 1)) {
logit_psi[1:n_sites, t + 1] <- gamma[t] + xmat %*% alphas #Get gamma colonization term as intercept in following years
}
# Ecological submodel: Define state conditional on parameters
for (i in 1:n_sites){
z[i,1] ~ dbern(psi[i,1]) #latent occupancy state in year 1 ~ bernoulli with prob based on our covariates and alpha01
psi[i,1] <- ilogit(logit_psi[i,1])
#Detection Model Year 1
for(j in 1:n_visits[i, 1]){
dets[i, j, 1] ~ dbern(z[i, 1] * p[i, j, 1])
p[i, j, 1] <- ilogit(beta0 +
beta1 * tmin[i, j, 1] +
beta2 * dayl[i, j, 1] +
beta3 * clut0[i, j, 1] +
beta4 * clut1[i, j, 1] +
beta5 * clut2[i, j, 1] +
beta6 * clut3[i, j, 1] +
beta7 * wind[i, j, 1])
}
#Occurrence Probs and Latent state for years 2 to n_years
for (t in 2:n_years){
z[i,t] ~ dbern(psi[i,t])
psi[i,t] <- ilogit(logit_psi[i, t] + z[i, (t-1)] * phi[t-1])
#Detection model for years 2 to n_years
for(j in 1:n_visits[i, t]){
dets[i, j, t] ~ dbern(z[i, t] * p[i, j, t])
p[i, j, t] <- ilogit(beta0 +
beta1 * tmin[i, j, t] +
beta2 * dayl[i, j, t] +
beta3 * clut0[i, j, t] +
beta4 * clut1[i, j, t] +
beta5 * clut2[i, j, t] +
beta6 * clut3[i, j, t] +
beta7 * wind[i, j, t])
}
}
}
# Derived parameters
# lambda measurements (trend)
lam.tot<-psi.hat[7]/psi.hat[1]
lam.avg<-mean(psi.hat[])
lam.tot.avg<-avg.psi[7]/avg.psi[1]
lam.avg.avg<-mean(avg.psi[])
#Original estimate for psi.hat (just the intercept of our original model)
psi.hat[1]<-ilogit(alpha01)
#Psi estimates for following years
for (t in 2:n_years){
col[t]<-ilogit(gamma[t-1])
surv[t]<-ilogit(gamma[t-1]+phi[t-1])
psi.hat[t]<-psi.hat[t-1]*surv[t]+(1-psi.hat[t-1])*col[t]
turnover[t]<-((1-psi.hat[t-1])*col[t])/((1-psi.hat[t-1])*col[t]+surv[t]*psi.hat[t-1])
}
#Avg psi at year 1
avg.psi[1] <- mean(psi[ ,1])
#Avg psi at year t
for(t in 2:n_years){
avg.psi[t] <- mean(psi[,t])
}
for (t in 1:n_years){
p.hat[t]<-ilogit(beta0)
}
}
